home *** CD-ROM | disk | FTP | other *** search
/ Amiga Plus Special 5 / Amiga Plus Sonderheft 1996 #5.iso / programme / povray / pov-ray_v2.2 / source / vect.c < prev    next >
C/C++ Source or Header  |  1995-04-29  |  25KB  |  994 lines
  1. /****************************************************************************
  2. *                vect.c
  3. *
  4. *  This file was written by Alexander Enzmann.  He wrote the code for
  5. *  4th-6th order shapes and generously provided us these enhancements.
  6. *
  7. *  from Persistence of Vision Raytracer
  8. *  Copyright 1993 Persistence of Vision Team
  9. *---------------------------------------------------------------------------
  10. *  NOTICE: This source code file is provided so that users may experiment
  11. *  with enhancements to POV-Ray and to port the software to platforms other 
  12. *  than those supported by the POV-Ray Team.  There are strict rules under
  13. *  which you are permitted to use this file.  The rules are in the file
  14. *  named POVLEGAL.DOC which should be distributed with this file. If 
  15. *  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  16. *  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
  17. *  Forum.  The latest version of POV-Ray may be found there as well.
  18. *
  19. * This program is based on the popular DKB raytracer version 2.12.
  20. * DKBTrace was originally written by David K. Buck.
  21. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  22. *
  23. *****************************************************************************/
  24.  
  25. #include "frame.h"
  26. #include "povproto.h"
  27. #include "vector.h"
  28.  
  29. #ifdef ABS
  30. #undef ABS
  31. #endif
  32.  
  33. #ifdef MAX
  34. #undef MAX
  35. #endif
  36.  
  37. extern int Shadow_Test_Flag;
  38. #undef EPSILON
  39. #define EPSILON 1.0e-10
  40. #define COEFF_LIMIT 1.0e-20
  41.  
  42. /*                  WARNING     WARNING    WARNING
  43.  
  44.    The following three constants have been defined so that quartic equations
  45.    will properly render on fpu/compiler combinations that only have 64 bit
  46.    IEEE precision.  Do not arbitrarily change any of these values.
  47.  
  48.    If you have a machine with a proper fpu/compiler combo - like a Mac with
  49.    a 68881, then use the native floating format (96 bits) and tune the
  50.    values for a little less tolerance: something like: factor1 = 1.0e15,
  51.    factor2 = -1.0e-7, factor3 = 1.0e-10.
  52.  
  53.    The meaning of the three constants are:
  54.       factor1 - the magnitude of difference between coefficients in a
  55.                 quartic equation at which the Sturmian root solver will
  56.         kick in.  The Sturm solver is quite a bit slower than
  57.         the algebraic solver, so the bigger this is, the faster
  58.         the root solving will go.  The algebraic solver is less
  59.         accurate so the bigger this is, the more likely you will
  60.         get bad roots.
  61.  
  62.       factor2 - Tolerance value that defines how close the quartic equation
  63.         is to being a square of a quadratic.  The closer this can
  64.         get to zero before roots disappear, the less the chance
  65.         you will get spurious roots.
  66.  
  67.       factor3 - Similar to factor2 at a later stage of the algebraic solver.
  68. */
  69. #define FUDGE_FACTOR1 1.0e12
  70. #define FUDGE_FACTOR2 -1.0e-5
  71. #define FUDGE_FACTOR3 1.0e-7
  72.  
  73. #define ABS(x) ((x) < 0.0 ? (0.0 - x) : (x))
  74. #define MAX(x,y) (x<y?y:x)
  75. #define TWO_PI 6.283185207179586476925286766560
  76. #define TWO_PI_3  2.0943951023931954923084
  77. #define TWO_PI_43 4.1887902047863909846168
  78. #define MAX_ITERATIONS 50
  79. #define MAXPOW 32
  80.  
  81.   typedef struct p {
  82.   int ord;
  83.   DBL coef[MAX_ORDER+1];
  84.   } 
  85. polynomial;
  86.  
  87. static int modp PARAMS((polynomial *u, polynomial *v, polynomial *r));
  88. int regula_falsa PARAMS((int order, DBL *coef, DBL a, DBL b, DBL *val));
  89. static int sbisect PARAMS((int np, polynomial *sseq, DBL min, DBL max, int atmin, int atmax, DBL *roots));
  90. int numchanges PARAMS((int np, polynomial *sseq, DBL a));
  91. DBL polyeval PARAMS((DBL x, int n, DBL *Coeffs));
  92. int buildsturm PARAMS((int ord, polynomial *sseq));
  93. int visible_roots PARAMS((int np, polynomial *sseq, int *atneg, int *atpos));
  94. static int difficult_coeffs PARAMS((int n, DBL *x));
  95.  
  96. /*
  97.  * modp
  98.  *
  99.  *   calculates the modulus of u(x) / v(x) leaving it in r, it
  100.  *  returns 0 if r(x) is a constant.
  101.  *  note: this function assumes the leading coefficient of v 
  102.  *   is 1 or -1
  103.  */
  104. static int modp(u, v, r)
  105. polynomial   *u, *v, *r;
  106.   {
  107.   int    i, k, j;
  108.   for (i=0;i<u->ord;i++)
  109.     r[i] = u[i];
  110.  
  111.   if (v->coef[v->ord] < 0.0) 
  112.     {
  113.     for (k = u->ord - v->ord - 1; k >= 0; k -= 2)
  114.       r->coef[k] = -r->coef[k];
  115.     for (k = u->ord - v->ord; k >= 0; k--)
  116.       for (j = v->ord + k - 1; j >= k; j--)
  117.       r->coef[j] = -r->coef[j] - r->coef[v->ord + k] * v->coef[j - k];
  118.     }
  119.   else 
  120.     {
  121.     for (k = u->ord - v->ord; k >= 0; k--)
  122.       for (j = v->ord + k - 1; j >= k; j--)
  123.       r->coef[j] -= r->coef[v->ord + k] * v->coef[j - k];
  124.     }
  125.  
  126.   k = v->ord - 1;
  127.   while (k >= 0 && fabs(r->coef[k]) < COEFF_LIMIT) 
  128.     {
  129.     r->coef[k] = 0.0;
  130.     k--;
  131.     }
  132.   r->ord = (k < 0) ? 0 : k;
  133.   return(r->ord);
  134.   }
  135.  
  136. /* Build the sturmian sequence for a polynomial */
  137. int buildsturm(ord, sseq)
  138. int      ord;
  139. polynomial   *sseq;
  140.   {
  141.   int i;
  142.   DBL f, *fp, *fc;
  143.   polynomial *sp;
  144.  
  145.   sseq[0].ord = ord;
  146.   sseq[1].ord = ord - 1;
  147.  
  148.   /* calculate the derivative and normalize the leading coefficient. */
  149.   f = fabs(sseq[0].coef[ord] * ord);
  150.   fp = sseq[1].coef;
  151.   fc = sseq[0].coef + 1;
  152.   for (i = 1; i <= ord; i++)
  153.     *fp++ = *fc++ * i / f;
  154.  
  155.   /* construct the rest of the Sturm sequence */
  156.   for (sp = sseq + 2;modp(sp - 2, sp - 1, sp); sp++) 
  157.     {
  158.     /* reverse the sign and normalize */
  159.     f = -fabs(sp->coef[sp->ord]);
  160.     for (fp = &sp->coef[sp->ord]; fp >= sp->coef; fp--)
  161.       *fp /= f;
  162.     }
  163.   sp->coef[0] = -sp->coef[0];   /* reverse the sign */
  164.   return(sp - sseq);
  165.   }
  166.  
  167. /* Find out how many visible intersections there are */
  168. int visible_roots(np, sseq, atzer, atpos)
  169. int np;
  170. polynomial *sseq;
  171. int *atzer, *atpos;
  172.   {
  173.   int atposinf, atzero;
  174.   polynomial *s;
  175.   DBL f, lf;
  176.  
  177.   atposinf = atzero = 0;
  178.   /* changes at positve infinity */
  179.   lf = sseq[0].coef[sseq[0].ord];
  180.   for (s = sseq + 1; s <= sseq + np; s++) 
  181.     {
  182.     f = s->coef[s->ord];
  183.     if (lf == 0.0 || lf * f < 0)
  184.       atposinf++;
  185.     lf = f;
  186.     }
  187.  
  188.   /* Changes at zero */
  189.   lf = sseq[0].coef[0];
  190.   for (s = sseq+1; s <= sseq + np; s++) 
  191.     {
  192.     f = s->coef[0];
  193.     if (lf == 0.0 || lf * f < 0)
  194.       atzero++;
  195.     lf = f;
  196.     }
  197.  
  198.   *atzer = atzero;
  199.   *atpos = atposinf;
  200.   return(atzero - atposinf);
  201.   }
  202.  
  203. /*
  204.  * numchanges
  205.  *
  206.  *   return the number of sign changes in the Sturm sequence in
  207.  * sseq at the value a.
  208.  */
  209. int numchanges(np, sseq, a)
  210. int      np;
  211. polynomial   *sseq;
  212. DBL   a;
  213.  
  214.   {
  215.   int      changes;
  216.   DBL   f, lf;
  217.   polynomial   *s;
  218.   changes = 0;
  219.   COOPERATE
  220.   lf = polyeval(a, sseq[0].ord, sseq[0].coef);
  221.   for (s = sseq + 1; s <= sseq + np; s++) 
  222.     {
  223.     f = polyeval(a, s->ord, s->coef);
  224.     if (lf == 0.0 || lf * f < 0)
  225.       changes++;
  226.     lf = f;
  227.     }
  228.   return(changes);
  229.   }
  230.  
  231. /*
  232.  * sbisect
  233.  *
  234.  *   uses a bisection based on the sturm sequence for the polynomial
  235.  * described in sseq to isolate intervals in which roots occur,
  236.  * the roots are returned in the roots array in order of magnitude.
  237.  
  238. Note: This routine has one severe bug: When the interval containing the
  239.       root [min, max] has a root at one of its endpoints, as well as one
  240.       within the interval, the root at the endpoint will be returned rather
  241.       than the one inside.
  242.  
  243.  */
  244. static int
  245. sbisect(np, sseq, min_value, max_value, atmin, atmax, roots)
  246. int np;
  247. polynomial *sseq;
  248. DBL min_value, max_value;
  249. int atmin, atmax;
  250. DBL *roots;
  251.   {
  252.   DBL  mid;
  253.   int  n1, n2, its, atmid;
  254.  
  255.   if ((atmin - atmax) == 1) 
  256.     {
  257.     /* first try using regula-falsa to find the root.  */
  258.     if (regula_falsa(sseq->ord, sseq->coef, min_value, max_value, roots))
  259.       return 1;
  260.     else 
  261.       {
  262.       /* That failed, so now find it by bisection */
  263.         for (its = 0; its < MAX_ITERATIONS; its++) 
  264.         {
  265.         mid = (min_value + max_value) / 2;
  266.         atmid = numchanges(np, sseq, mid);
  267.         if (fabs(mid) > EPSILON) 
  268.           {
  269.           if (fabs((max_value - min_value) / mid) < EPSILON) 
  270.             {
  271.             roots[0] = mid;
  272.             return 1;
  273.             }
  274.           }
  275.         else if (fabs(max_value - min_value) < EPSILON) 
  276.           {
  277.           roots[0] = mid;
  278.           return 1;
  279.           }
  280.         else if ((atmin - atmid) == 0)
  281.           min_value = mid;
  282.         else
  283.           max_value = mid;
  284.         }
  285.       /* Bisection took too long - just return what we got */
  286.       roots[0] = mid;
  287.       return 1;
  288.       }
  289.     }
  290.  
  291.   /* There is more than one root in the interval.
  292.       Bisect to find new intervals */
  293.   for (its = 0; its < MAX_ITERATIONS; its++) 
  294.     {
  295.     mid = (min_value + max_value) / 2;
  296.     atmid = numchanges(np, sseq, mid);
  297.     n1 = atmin - atmid;
  298.     n2 = atmid - atmax;
  299.     if (n1 != 0 && n2 != 0) 
  300.       {
  301.       n1 = sbisect(np, sseq, min_value, mid, atmin, atmid, roots);
  302.       n2 = sbisect(np, sseq, mid, max_value, atmid, atmax, &roots[n1]);
  303.       return n1 + n2;
  304.       }
  305.     if (n1 == 0)
  306.       min_value = mid;
  307.     else
  308.       max_value = mid;
  309.     }
  310.  
  311.   /* Took too long to bisect.  Just return what we got. */
  312.   roots[0] = mid;
  313.   return 1;
  314.   }
  315.  
  316. DBL polyeval(x, n, Coeffs)
  317. DBL x, *Coeffs;
  318. int n;
  319.   {
  320.   register int i;
  321.   DBL val;
  322.   val = Coeffs[n];
  323.   for (i=n-1;i>=0;i--) val = val * x + Coeffs[i];
  324.   return val;
  325.   }
  326.  
  327. /* Close in on a root by using regula-falsa */
  328. int regula_falsa(order, coef, a, b, val)
  329. int order;
  330. DBL *coef;
  331. DBL a, b, *val;
  332.   {
  333.   int its;
  334.   DBL fa, fb, x, fx, lfx;
  335.  
  336.   fa = polyeval(a, order, coef);
  337.   fb = polyeval(b, order, coef);
  338.  
  339.   if (fa * fb > 0.0)
  340.     return 0;
  341.  
  342.   if (fabs(fa) < COEFF_LIMIT) 
  343.     {
  344.     *val = a;
  345.     return 1;
  346.     }
  347.  
  348.   if (fabs(fb) < COEFF_LIMIT) 
  349.     {
  350.     *val = b;
  351.     return 1;
  352.     }
  353.  
  354.   lfx = fa;
  355.  
  356.   COOPERATE
  357.   for (its = 0; its < MAX_ITERATIONS; its++) 
  358.     {
  359.     x = (fb * a - fa * b) / (fb - fa);
  360.     fx = polyeval(x, order, coef);
  361.  
  362.     if (fabs(x) > EPSILON) 
  363.       {
  364.       if (fabs(fx / x) < EPSILON) 
  365.         {
  366.         *val = x;
  367.         return 1;
  368.         }
  369.       }
  370.     else if (fabs(fx) < EPSILON) 
  371.       {
  372.       *val = x;
  373.       return 1;
  374.       }
  375.  
  376.     if (fa < 0)
  377.       if (fx < 0) 
  378.       {
  379.       a = x;
  380.       fa = fx;
  381.       if ((lfx * fx) > 0)
  382.         fb /= 2;
  383.       }
  384.       else 
  385.       {
  386.       b = x;
  387.       fb = fx;
  388.       if ((lfx * fx) > 0)
  389.         fa /= 2;
  390.       }
  391.     else if (fx < 0) 
  392.       {
  393.       b = x;
  394.       fb = fx;
  395.       if ((lfx * fx) > 0)
  396.         fa /= 2;
  397.       }
  398.     else 
  399.       {
  400.       a = x;
  401.       fa = fx;
  402.       if ((lfx * fx) > 0)
  403.         fb /= 2;
  404.       }
  405.     if (fabs(b-a) < EPSILON) 
  406.       {
  407.       /* Check for underflow in the domain */
  408.       *val = x;
  409.       return 1;
  410.       }
  411.     lfx = fx;
  412.     }
  413.   return 0;
  414.   }
  415.  
  416. /*
  417.    Solve the quadratic equation:
  418.               x[0] * x^2 + x[1] * x + x[2] = 0.
  419.  
  420.    The value returned by this function is the number of real roots.
  421.    The roots themselves are returned in y[0], y[1].
  422. */
  423. int solve_quadratic(x, y)
  424. DBL *x, *y;
  425.   {
  426.   DBL d, t, a, b, c;
  427.   a = x[0];
  428.   b = -x[1];
  429.   c = x[2];
  430.   if (a == 0.0) 
  431.     {
  432.     if (b == 0.0)
  433.       return 0;
  434.     y[0] = c / b;
  435.     return 1;
  436.     }
  437.   d = b * b - 4.0 * a * c;
  438.   if (d < 0.0)
  439.     return 0;
  440.   else if (fabs(d) < COEFF_LIMIT) 
  441.     {
  442.     y[0] = 0.5 * b / a;
  443.     return 1;
  444.     }
  445.   d = sqrt(d);
  446.   t = 2.0 * a;
  447.   y[0] = (b + d) / t;
  448.   y[1] = (b - d) / t;
  449.   return 2;
  450.   }
  451.  
  452. /*
  453.    Solve the cubic equation:
  454.  
  455.       x[0] * x^3 + x[1] * x^2 + x[2] * x + x[3] = 0.
  456.  
  457.    The result of this function is an integer that tells how many real
  458.    roots exist.  Determination of how many are distinct is up to the
  459.    process that calls this routine.  The roots that exist are stored
  460.    in (y[0], y[1], y[2]).
  461.  
  462.    Note: this function relies very heavily on trigonometric functions and
  463.    the square root function.  If an alternative solution is found that does
  464.    not rely on transcendentals this code will be replaced.
  465. */
  466. int solve_cubic(x, y)
  467. DBL *x, *y;
  468.   {
  469.   DBL Q, R, Q3, R2, sQ, d, an, theta;
  470.   DBL A2, a0, a1, a2, a3;
  471.   a0 = x[0];
  472.   if (a0 == 0.0) return solve_quadratic(&x[1], y);
  473.   else if (a0 != 1.0) 
  474.     {
  475.     a1 = x[1] / a0;
  476.     a2 = x[2] / a0;
  477.     a3 = x[3] / a0;
  478.     }
  479.   else 
  480.     {
  481.     a1 = x[1];
  482.     a2 = x[2];
  483.     a3 = x[3];
  484.     }
  485.   A2 = a1 * a1;
  486.   Q = (A2 - 3.0 * a2) / 9.0;
  487.   R = (2.0 * A2 * a1 - 9.0 * a1 * a2 + 27.0 * a3) / 54.0;
  488.   Q3 = Q * Q * Q;
  489.   R2 = R * R;
  490.   d = Q3 - R2;
  491.   an = a1 / 3.0;
  492.   if (d >= 0.0) 
  493.     {
  494.     /* Three real roots. */
  495.     d = R / sqrt(Q3);
  496.     theta = acos(d) / 3.0;
  497.     sQ = -2.0 * sqrt(Q);
  498.     y[0] = sQ * cos(theta) - an;
  499.     y[1] = sQ * cos(theta + TWO_PI_3) - an;
  500.     y[2] = sQ * cos(theta + TWO_PI_43) - an;
  501.     return 3;
  502.     }
  503.   else 
  504.     {
  505.     sQ = pow(sqrt(R2 - Q3) + ABS(R), 1.0 / 3.0);
  506.     if (R < 0)
  507.       y[0] = (sQ + Q / sQ) - an;
  508.     else
  509.       y[0] = -(sQ + Q / sQ) - an;
  510.     return 1;
  511.     }
  512.   }
  513.  
  514. /* Test to see if any coeffs are more than 6 orders of magnitude
  515.    larger than the smallest */
  516. static int
  517. difficult_coeffs(n, x)
  518. int n;
  519. DBL *x;
  520.   {
  521.   int i, flag = 0;
  522.   DBL t, biggest;
  523.  
  524.   biggest = fabs(x[0]);
  525.   for (i=1;i<=n;i++) 
  526.     {
  527.     t = fabs(x[i]);
  528.     if (t > biggest)
  529.       biggest = t;
  530.     }
  531.  
  532.   /* Everything is zero no sense in doing any more */
  533.   if (biggest == 0.0) return 0;
  534.  
  535.   for (i=0;i<=n;i++)
  536.     if (x[i] != 0.0)
  537.       if (fabs(biggest / x[i]) > FUDGE_FACTOR1) 
  538.       {
  539.       x[i] = 0.0;
  540.       flag = 1;
  541.       }
  542.  
  543.   return flag;
  544.   }
  545.  
  546. #ifdef TEST_SOLVER
  547. /* The old way of solving quartics algebraically */
  548. /* This is an adaptation of the method of Lodovico Ferrari (Circa 1731). */
  549. int solve_quartic(x, results)
  550. DBL *x, *results;
  551.   {
  552.   DBL cubic[4], roots[3];
  553.   DBL a0, a1, y, d1, x1, t1, t2;
  554.   DBL c0, c1, c2, c3, c4, d2, q1, q2;
  555.   int i;
  556.  
  557.   /* Figure out the size difference between coefficients */
  558.   if (difficult_coeffs(4, x)) 
  559.     {
  560.     if (fabs(x[0]) < COEFF_LIMIT)
  561.       if (fabs(x[1]) < COEFF_LIMIT)
  562.         return solve_quadratic(&x[2], results);
  563.       else
  564.         return solve_cubic(&x[1], results);
  565.     else
  566.       return polysolve(4, x, results);
  567.     }
  568.  
  569.   c0 = x[0];
  570.   if (fabs(c0) < COEFF_LIMIT)
  571.     return solve_cubic(&x[1], results);
  572.   else if (fabs(x[4]) < COEFF_LIMIT) 
  573.     {
  574.     return solve_cubic(x, results);
  575.     }
  576.   else if (c0 != 1.0) 
  577.     {
  578.     c1 = x[1] / c0;
  579.     c2 = x[2] / c0;
  580.     c3 = x[3] / c0;
  581.     c4 = x[4] / c0;
  582.     }
  583.   else 
  584.     {
  585.     c1 = x[1];
  586.     c2 = x[2];
  587.     c3 = x[3];
  588.     c4 = x[4];
  589.     }
  590.  
  591.     /* The first step is to take the original equation:
  592.  
  593.      x^4 + b*x^3 + c*x^2 + d*x + e = 0
  594.  
  595.       and rewrite it as:
  596.  
  597.      x^4 + b*x^3 = -c*x^2 - d*x - e,
  598.  
  599.       adding (b*x/2)^2 + (x^2 + b*x/2)y + y^2/4 to each side gives a
  600.       perfect square on the lhs:
  601.  
  602.      (x^2 + b*x/2 + y/2)^2 = (b^2/4 - c + y)x^2 + (b*y/2 - d)x + y^2/4 - e
  603.  
  604.       By choosing the appropriate value for y, the rhs can be made a perfect
  605.       square also.  This value is found when the rhs is treated as a quadratic
  606.       in x with the discriminant equal to 0.  This will be true when:
  607.  
  608.      (b*y/2 - d)^2 - 4.0 * (b^2/4 - c*y)*(y^2/4 - e) = 0, or
  609.  
  610.      y^3 - c*y^2 + (b*d - 4*e)*y - b^2*e + 4*c*e - d^2 = 0.
  611.  
  612.       This is called the resolvent of the quartic equation.  */
  613.  
  614.     a0 = 4.0 * c4;
  615.   cubic[0] = 1.0;
  616.   cubic[1] = -1.0 * c2;
  617.   cubic[2] = c1 * c3 - a0;
  618.   cubic[3] = a0 * c2 - c1 * c1 * c4 - c3 * c3;
  619.   i = solve_cubic(&cubic[0], &roots[0]);
  620.   if (i > 0)
  621.     y = roots[0];
  622.   else
  623.     return 0;
  624.  
  625.   /* What we are left with is a quadratic squared on the lhs and a
  626.       linear term on the right.  The linear term has one of two signs,
  627.       take each and add it to the lhs.  The form of the quartic is now:
  628.  
  629.      a' = b^2/4 - c + y,    b' = b*y/2 - d, (from rhs quadritic above)
  630.  
  631.      (x^2 + b*x/2 + y/2) = +sqrt(a'*(x + 1/2 * b'/a')^2), and
  632.      (x^2 + b*x/2 + y/2) = -sqrt(a'*(x + 1/2 * b'/a')^2).
  633.  
  634.       By taking the linear term from each of the right hand sides and
  635.       adding to the appropriate part of the left hand side, two quadratic
  636.       formulas are created.  By solving each of these the four roots of
  637.       the quartic are determined.
  638.    */
  639.   i = 0;
  640.   a0 = c1 / 2.0;
  641.   a1 = y / 2.0;
  642.  
  643.   t1 = a0 * a0 - c2 + y;
  644.   if (t1 < 0.0) 
  645.     {
  646.     if (t1 > FUDGE_FACTOR2)
  647.       t1 = 0.0;
  648.     else 
  649.       {
  650.       /* First Special case, a' < 0 means all roots are complex. */
  651.         return 0;
  652.       }
  653.     }
  654.   if (t1 < FUDGE_FACTOR3) 
  655.     {
  656.     /* Second special case, the "x" term on the right hand side above
  657.      has vanished.  In this case:
  658.         (x^2 + b*x/2 + y/2) = +sqrt(y^2/4 - e), and
  659.         (x^2 + b*x/2 + y/2) = -sqrt(y^2/4 - e).  */
  660.     t2 = a1 * a1 - c4;
  661.     if (t2 < 0.0) 
  662.       {
  663.       return 0;
  664.       }
  665.     x1 = 0.0;
  666.     d1 = sqrt(t2);
  667.     }
  668.   else 
  669.     {
  670.     x1 = sqrt(t1);
  671.     d1 = 0.5 * (a0 * y - c3) / x1;
  672.     }
  673.   /* Solve the first quadratic */
  674.     q1 = -a0 - x1;
  675.   q2 = a1 + d1;
  676.   d2 = q1 * q1 - 4.0 * q2;
  677.   if (d2 >= 0.0) 
  678.     {
  679.     d2 = sqrt(d2);
  680.     results[0] = 0.5 * (q1 + d2);
  681.     results[1] = 0.5 * (q1 - d2);
  682.     i = 2;
  683.     }
  684.   /* Solve the second quadratic */
  685.   q1 = q1 + x1 + x1;
  686.   q2 = a1 - d1;
  687.   d2 = q1 * q1 - 4.0 * q2;
  688.   if (d2 >= 0.0) 
  689.     {
  690.     d2 = sqrt(d2);
  691.     results[i++] = 0.5 * (q1 + d2);
  692.     results[i++] = 0.5 * (q1 - d2);
  693.     }
  694.   return i;
  695.   }
  696. #else
  697.   /* Solve a quartic using the method of Francois Vieta (Circa 1735) */
  698.   int
  699.   solve_quartic(x, results)
  700.     DBL *x, *results;
  701.     {
  702.     DBL cubic[4], roots[3];
  703.     DBL c12, z, p, q, q1, q2, r, d1, d2;
  704.     DBL c0, c1, c2, c3, c4;
  705.     int i;
  706.  
  707.     /* Figure out the size difference between coefficients */
  708.     if (difficult_coeffs(4, x)) 
  709.       {
  710.       if (fabs(x[0]) < COEFF_LIMIT)
  711.         if (fabs(x[1]) < COEFF_LIMIT)
  712.           return solve_quadratic(&x[2], results);
  713.         else
  714.           return solve_cubic(&x[1], results);
  715.       else
  716.         return polysolve(4, x, results);
  717.       }
  718.  
  719.     /* See if the high order term has vanished */
  720.     c0 = x[0];
  721.     if (fabs(c0) < COEFF_LIMIT)
  722.       return solve_cubic(&x[1], results);
  723.  
  724.     /* See if the constant term has vanished */
  725.     if (fabs(x[4]) < COEFF_LIMIT)
  726.       return solve_cubic(x, results);
  727.  
  728.     /* Make sure the quartic has a leading coefficient of 1.0 */
  729.     if (c0 != 1.0) 
  730.       {
  731.       c1 = x[1] / c0;
  732.       c2 = x[2] / c0;
  733.       c3 = x[3] / c0;
  734.       c4 = x[4] / c0;
  735.       }
  736.     else 
  737.       {
  738.       c1 = x[1];
  739.       c2 = x[2];
  740.       c3 = x[3];
  741.       c4 = x[4];
  742.       }
  743.  
  744.       /* Compute the cubic resolvant */
  745.       c12 = c1 * c1;
  746.     p = -0.375 * c12 + c2;
  747.     q = 0.125 * c12 * c1 - 0.5 * c1 * c2 + c3;
  748.     r = -0.01171875 * c12 * c12 + 0.0625 * c12 * c2 - 0.25 * c1 * c3 + c4;
  749.  
  750.     cubic[0] = 1.0;
  751.     cubic[1] = -0.5 * p;
  752.     cubic[2] = -r;
  753.     cubic[3] = 0.5 * r * p - 0.125 * q * q;
  754.     i = solve_cubic(cubic, roots);
  755.     if (i > 0)
  756.       z = roots[0];
  757.     else
  758.       return 0;
  759.  
  760.     d1 = 2.0 * z - p;
  761.  
  762.     if (d1 < 0.0) 
  763.       {
  764.       if (d1 > -EPSILON)
  765.         d1 = 0.0;
  766.       else
  767.         return 0;
  768.       }
  769.     if (d1 < EPSILON) 
  770.       {
  771.       d2 = z * z - r;
  772.       if (d2 < 0.0)
  773.         return 0;
  774.       d2 = sqrt(d2);
  775.       }
  776.     else 
  777.       {
  778.       d1 = sqrt(d1);
  779.       d2 = 0.5 * q / d1;
  780.       }
  781.  
  782.       /* Set up useful values for the quadratic factors */
  783.       q1 = d1 * d1;
  784.     q2 = -0.25 * c1;
  785.     i = 0;
  786.  
  787.     /* Solve the first quadratic */
  788.     p = q1 - 4.0 * (z - d2);
  789.     if (p == 0)
  790.       results[i++] = -0.5 * d1 - q2;
  791.     else if (p > 0) 
  792.       {
  793.       p = sqrt(p);
  794.       results[i++] = -0.5 * (d1 + p) + q2;
  795.       results[i++] = -0.5 * (d1 - p) + q2;
  796.       }
  797.     /* Solve the second quadratic */
  798.     p = q1 - 4.0 * (z + d2);
  799.     if (p == 0)
  800.       results[i++] = 0.5 * d1 - q2;
  801.     else if (p > 0) 
  802.       {
  803.       p = sqrt(p);
  804.       results[i++] = 0.5 * (d1 + p) + q2;
  805.       results[i++] = 0.5 * (d1 - p) + q2;
  806.       }
  807.     return i;
  808.     }
  809. #endif
  810.  
  811. /* Root solver based on the Sturm sequences for a polynomial. */
  812. int polysolve(order, Coeffs, roots)
  813. int order;
  814. DBL *Coeffs, *roots;
  815.   {
  816.   polynomial sseq[MAX_ORDER+1];
  817.   DBL min_value, max_value;
  818.   int i, nroots, np, atmin, atmax;
  819.  
  820.   /* Put the coefficients into the top of the stack. */
  821.   for (i=0;i<=order;i++)
  822.     sseq[0].coef[order-i] = Coeffs[i];
  823.  
  824.   /* Build the Sturm sequence */
  825.   np = buildsturm(order, &sseq[0]);
  826.  
  827.   /* Get the total number of visible roots */
  828.   if ((nroots = visible_roots(np, sseq, &atmin, &atmax)) == 0)
  829.     return 0;
  830.  
  831.   /* Bracket the roots */
  832.   min_value = 0.0;
  833.   max_value = Max_Distance;
  834.  
  835.   atmin = numchanges(np, sseq, min_value);
  836.   atmax = numchanges(np, sseq, max_value);
  837.   nroots = atmin - atmax;
  838.   if (nroots == 0) return 0;
  839.  
  840.   /* perform the bisection. */
  841.   return sbisect(np, sseq, min_value, max_value, atmin, atmax, roots);
  842.   }
  843.  
  844. #ifdef TEST     /* This is not used anywhere or tested.  Interesting? */
  845.  
  846. #define MAX_POLYGON_SIDES 8
  847. #define Crossing_Point(x1, y1, x2, y2) (x1 - y1 * (x2 - x1) / (y2 - y1))
  848.  
  849. /* See if "Ray" intersects the polygon defined by the coordinate list "vertices". */
  850. int Intersect_Polygon(Ray, vertex_count, vertices, n, d, Depth, Intersect_Point)
  851. RAY *Ray;
  852. int vertex_count;
  853. VECTOR *vertices, *n, *Intersect_Point;
  854. DBL d, *Depth;
  855.   {
  856.   DBL s, t, x, y, z;
  857.   int Sign_Holder, Next_Sign, Crossings;
  858.   int i, this_vertex, next_vertex;
  859.  
  860.   static DBL temp_x[MAX_POLYGON_SIDES],
  861.   temp_y[MAX_POLYGON_SIDES];
  862.  
  863.   /* Calculate the point of intersection and the depth. */
  864.   VDot(s, Ray->Direction, *n);
  865.   if (s == 0.0)
  866.     return 0;
  867.   VDot(t, Ray->Initial, *n);
  868.   *Depth = 0.0 - (d + t) / s;
  869.   if (*Depth < Small_Tolerance)
  870.     return 0;
  871.   VScale(*Intersect_Point, Ray->Direction, *Depth);
  872.   VAdd(*Intersect_Point, *Intersect_Point, Ray->Initial);
  873.  
  874.   /* Map the intersection point and the triangle onto a plane. */
  875.   x = ABS(n->x); y = ABS(n->y); z = ABS(n->z);
  876.   if (x>y)
  877.     if (x>z)
  878.       for (i=0;i<vertex_count;i++) 
  879.     {
  880.     temp_x[i] = vertices[i].y - Intersect_Point->y;
  881.     temp_y[i] = vertices[i].z - Intersect_Point->z;
  882.     }
  883.   else
  884.     for (i=0;i<vertex_count;i++) 
  885.     {
  886.     temp_x[i] = vertices[i].x - Intersect_Point->x;
  887.     temp_y[i] = vertices[i].y - Intersect_Point->y;
  888.     }
  889.   else if (y>z)
  890.     for (i=0;i<vertex_count;i++) 
  891.     {
  892.     temp_x[i] = vertices[i].x - Intersect_Point->x;
  893.     temp_y[i] = vertices[i].z - Intersect_Point->z;
  894.     }
  895.   else
  896.     for (i=0;i<vertex_count;i++) 
  897.     {
  898.     temp_x[i] = vertices[i].x - Intersect_Point->x;
  899.     temp_y[i] = vertices[i].y - Intersect_Point->y;
  900.     }
  901.  
  902.   /* Determine crossing count */
  903.   Crossings = 0;
  904.   if (temp_y[0] < 0) Sign_Holder = -1;
  905.   else Sign_Holder = 1;
  906.  
  907.     for (i=0;i<vertex_count;i++) 
  908.     {
  909.     /* Start of winding tests, test the segment from v1 to v2 */
  910.     this_vertex = i;
  911.     next_vertex = (i + 1) % vertex_count;
  912.     if (temp_y[next_vertex] < 0) Next_Sign = -1;
  913.     else Next_Sign = 1;
  914.     if (Sign_Holder != Next_Sign) 
  915.       {
  916.       if (temp_x[this_vertex] > 0.0) 
  917.         {
  918.         if (temp_x[next_vertex] > 0.0)
  919.           Crossings++;
  920.         else if (Crossing_Point(temp_x[this_vertex], temp_y[this_vertex],
  921.           temp_x[next_vertex], temp_y[next_vertex]) > 0.0)
  922.           Crossings++;
  923.         }
  924.       else if (temp_x[next_vertex] > 0.0) 
  925.         {
  926.         if (Crossing_Point(temp_x[this_vertex], temp_y[this_vertex],
  927.           temp_x[next_vertex], temp_y[next_vertex]) > 0.0)
  928.           Crossings++;
  929.         }
  930.       Sign_Holder = Next_Sign;
  931.       }
  932.     }
  933.  
  934.   return (Crossings&1); /* Inside only if # of crossings odd. */
  935.   }
  936.  
  937. /* Translate screen coordinates into world coordinates. */
  938. void World_Coordinate (Viewpoint, Ray, width, height, x, y)
  939. VIEWPOINT *Viewpoint;
  940. RAY *Ray;
  941. int width, height;
  942. DBL x, y;
  943.   {
  944.   DBL scalex, scaley;
  945.   VECTOR V1, V2;
  946.  
  947.   /* Convert the X Coordinate to be a DBL from -0.5 to 0.5 */
  948.   scalex = (x - (DBL)width / 2.0) / (DBL)width;
  949.   /* Convert the Y Coordinate to be a DBL from -0.5 to 0.5 */
  950.   scaley = (((DBL)(height - 1) - y) - (DBL)height / 2.0) / (DBL)height;
  951.   /* Compute the direction of the screen point from the viewpoint */
  952.   VScale (V1, Viewpoint->Up, scaley);
  953.   VScale (V2, Viewpoint->Right, scalex);
  954.   VAdd (Ray->Direction, V1, V2);
  955.   VAdd (Ray->Direction, Ray->Direction, Viewpoint->Direction);
  956.   VNormalize (Ray->Direction, Ray->Direction);
  957.   }
  958. /* Uncomment these two declarations if your compiler needs them */
  959. /* They give Microsoft C an out of macro expansion space error */
  960. /* void show_univariate_poly PARAMS((char *label, int order,DBL *Coeffs));  */
  961. /* void show_points PARAMS((char *label,int count,DBL *point_list);  */
  962.  
  963. void show_univariate_poly(label, order, Coeffs)
  964. char *label;
  965. int order;
  966. DBL *Coeffs;
  967.   {
  968.   int i;
  969.   printf("%s", label);
  970.   for (i=0;i<=order;i++) 
  971.     {
  972.     printf("%.2lf x^%d", Coeffs[i], order-i);
  973.     if (i==order) printf("\n");
  974.     else printf(" + ");
  975.     }
  976.   }
  977.  
  978.   void show_points(label, count, point_list)
  979.     char *label;
  980. int count;
  981. DBL *point_list;
  982.   {
  983.   int i;
  984.   printf("%s", label);
  985.   for (i=0;i<count;i++) 
  986.     {
  987.     printf("%lf", point_list[i]);
  988.     if (i==(count-1)) printf("\n");
  989.     else printf(", ");
  990.     }
  991.   }
  992.  
  993. #endif
  994.